Expansion dynamics of Lennard-Jones systems 
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The dynamics of the expansion of a Lennard-Jones system, initially confined at high density and 
subsequently expanding freely in the vacuum, is confronted to an expanding statistical ensemble, 
derived in the diluted quasi-ideal Boltzmann approximation. The description proves to be fairly 
accurate at predicting average one-body global observables, but important deviations are observed 
in the configuration-space structure of the events. Possible implications for finite expanding physical 
systems are outlined. 
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INTRODUCTION 

When a finite isolated unbound system decays in the 
vacuum, its decay pattern is often characterized by an or- 
dered kinetic component: the collective flow. This is the 
case in clusters dissociation induced by photo-ionization 
0) B B 0] or charge transfer collisions 0, Here, the 
condensed matter bulk limit encourages, in principle, an 
interpretation in terms of liquid-vapor transition. How- 
ever, the obvious fact that all vapor flows out and no 
vapor comes back in, makes it difficult to push the anal- 
ogy further, and the dynamical evaporation can only be 
interpreted thermodynamically [l], [2, H| by making use of 
a time-dependent temperature within the concept of an 
evaporative ensemble 0, Q • Flow is also a basic feature 
of heavy ion collisions, where the products of fragmenta- 
tion reactions show a velocity preferentially oriented in 
the radial direction [9(. If in the Fermi energy regime 
and in the associated multi-fragmentation phase transi- 
tion these collective flows may be only a perturbation 
in the global energetics, this is not true at SIS ener- 
gies (between 0.2 and 2.0 GeV/nucleon), where they are 
likely to influence light cluster formation by coalescence 
[10]. In the ultra-relativistic regime the ordered and dis- 
ordered motions become comparable in magnitude [ill ], 
and collective flows are believed to play an essential role 
in the characteristics of the transition to the quark-gluon 
plasma observed in the Relativistic Heavy Ion Collider 

data [mmm. 

In the initial stage of a nuclear collision the complex- 
ity of the dynamics is such that a statistical analysis of 
the system might prove useful, even at energies as low 
as in the Fermi regime. If this equilibrium stage occurs 
at high density [14j, the final-state interaction may still 
be important in the subsequent evolution, and it is not 



clear how the final partitions at the freeze-out stage will 
be modified. In a previous paper [lH ] we have addressed 
this issue at a classical level through molecular dynam- 
ics simulations of a schematic Lennard-Jones system [l6[ , 
initially thermalized at high (supercritical) density, and 
subsequently freely expanding in the vacuum. We have 
shown that the dynamics of the expansion leads to a con- 
siderable increase of fluctuations. By the time when par- 
titions are settled and the formed pre-fragments cease to 
interact (freeze-out), these fluctuations are qualitatively 
similar to the ones expected from a thermal system at re- 
duced subcritical density. A naturally arising question is 
then whether the expanded system at freeze-out can still 
be treated as a statistical equilibrium at a lower density. 

Our numerical experiments can only address in a very 
partial and incomplete way this important question, 
which does probably not have a unique answer. In the 
evaporation regime accessed in cluster experiments the 
above-mentioned freeze-out configuration does not exist, 
and a statistical treatment equivalent to the dynamical 
process has to be adapted to the time window of the 
experiment [l[. Conversely in a heavy ion collision, time 
scales are so short that well defined freeze-out time(s) can 
be identified. However, such collisions may correspond to 
a dynamical process very different from the free expan- 
sion of a dense system of classical particles. In particu- 
lar, it is important to stress that macroscopic statistical 
models in all energy regimes 17, [l8| suppose that the 
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statistical hypothesis applies at the freeze-out time, and 
in this sense they by-pass by construction the problem of 
the out-of-equilibrium evolution up to freeze-out. 

In different physical situations flow appears to settle 
early in the dynamics. Notably, this is the case of cen- 
tral nuclear collisions, where flow is associated to an ini- 
tial compression of the dinuclcar system. In such a sit- 
uation equilibrium, if ever reached, is approached when 
the system is still strongly interacting, and we may ex- 
pect the successive dynamics to deeply modify the system 
configurations. In these cases, the model of an initially 
equilibrated dense molecular system, which then is freely 
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expanding in the vacuum, may bear some pertinent in- 
formation. 

In this paper we compare freeze-out configurations of 
a freely expanding Lennard- Jones system with different 
equilibrium models for the same system. We show that, if 
statistical models are reasonably correct as far as average 
quantities are concerned, important differences can be 
seen in the fragmentation patterns. In particular, close 
to the liquid-gas phase transition unstable configurations, 
inaccessible to equilibrium models, appear to dominate 
the dynamics of the expansion. 



new velocity drawn from a Maxwell-Boltzmann distribu- 
tion at the desired temperature /3" 1 . Between stochastic 
collisions, the system evolves at constant energy. It has 
been proved that, under some general conditions (l9l |. 
the constant energy shells are visited according to their 
Boltzmann weights, which in turn implies that the en- 
semble of configurations at different times constitutes a 
canonical ensemble at the thermostat temperature, dis- 
tributed as 

„ M . _Lexp-^[<aJ- + ^1">)> 
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THE ISOBAR MICROCANONICAL ENSEMBLE 

The system under study is composed of excited drops 
made up of N = 147 particles interacting via a Lennard- 
Jones (LJ) 6—12 potential VLj(r) with a cutoff radius 
r c = 3o\ Energies are measured in units of the potential 
well (e) , a characterizes the radius of a particle and m is 
its mass. We adopt adimensional units for energy, length, 
and time such that e = a = 1, to = y/ er 2 m/48e. The 
initial condition is given by the microcanonical isobar 
statistical ensemble described by the probability for each 
microstate (n) 



N 



(4) 
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Microcanonical ensembles are extracted at different en- 
ergies by sorting the events of the canonical distributions 
according to their energy, excluding the contribution of 
the confining potential. These dense configurations cor- 
respond for all energies to the supercritical part of the 
Lennard- Jones phase diagram (20j . 



COMPARISON WITH THE FREE EXPANDING 
SYSTEM 



(n) 
Po 



where r| n ' is the position of particle i within the mi- 
crostate (n), H^j is the corresponding Lennard- Jones 
energy, A is a Lagrange multiplier constraining a finite 
size, and 
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W X (E) = 5>p[-A5>^) 2 ]<5 (e - H™) 

(n) i=l 



(2) 



is the associated density of states or partition sum. The 
distribution eq. ([I]) is the minimum biased probability dis- 
tribution for an isolated finite system with a finite size 
measured by its mean square radius [2~fl | 



N 

t=l 



(3) 



where the average is taken over microstates. To generate 
the statistical ensemble eq.fl]), we numerically proceed 
as follows. A harmonic potential with spring constant 
k = moj 2 , with lj — O.Uq 1 is added to the Hamiltonian, 
and the system is coupled to a thermostat using the An- 
dersen technique (l9j to achieve equilibrium inside the os- 
cillator. In brief, this is attained by stochastic impulsive 
forces that act occasionally on randomly selected parti- 
cles. After each collision the particle is endowed with a 



Choice of the comparison time 

The equilibrated configurations are let evolve in the 
vacuum for a time long enough that the chemical com- 
position of the system is settled (freeze-out time) [27J • An 
interesting observable to study the freeze-out properties 
of the system is given by the normalized kinetic energy 
fluctuation Ak — N o\ / (K) 2 . This quantity shows a 
saturation when the sharing between kinetic and poten- 
tial energy finishes, which physically corresponds to the 
formation of surfaces within the system, which in turn 
implies a chemical stabilization of the cluster properties 
|15j . The behavior of this fluctuation with time is dis- 
played for four representative energy states in the upper 
part of Figure 1 . 

During the time evolution, the initially disordered ki- 
netic energy is partially converted into collective motion, 
defined as 
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The internal energy at each time is correspondingly given 
by E int = K int + V LJ , where 



N 

E 

i=l 



2m ( 
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Pi 
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part of the Lennard-Jones energy. 



and = Ei 
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is the interaction 
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The time evolution of the system with different initial 
energies is represented in the (E int ) versus (R 2 ) plane 
in the lower part of Figure 1. We can see that the time 
evolution corresponds to a rapid expansion and a devel- 
opment of collective flow which enlarges with increas- 
ing total energy. If at the entropy-saturation time or 
freeze-out time t / the system were still close to a statis- 
tical equilibrium, it should be described by eq.fll]) with 
E = (Ei nt )(tf) and a new volume constraint uj such that 
the average over the statistical ensemble (JXJ) of the mean 
square radius is (R 2 ) = (R 2 )(tf). The mean square ra- 
dius of an equilibrated system with the looser constraint 
uj = O.Ol^ 1 is also represented in the lower part of Figure 
1. The points where this curve crosses the time evolution 
of the freely expanding systems give the times we have 
chosen to compare the statistical and dynamical ensem- 
bles. As it can be observed from the upper part of the 
figure, for all chosen energies these times are long enough 
for the configurations to be safely considered as frozen. 
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FIG. 1: (Color online) Lower part: mean square radius 
as a function of the internal energy of a confined LJ sys- 
tem in a harmonic trap with uj — O.Olt^ 1 (full line) and 
uj — Q.Uq (dashed line), compared to LJ systems in free 
expansion, at different total energies and times between 0to 
and tf = 250to (symbols). Upper part: kinetic energy fluctua- 
tion Ak (see text) as a function of time for the free expanding 
systems at 4 different energies. The large circles in both pic- 
tures correspond to the time at which the confined and the 
free expanding system have the same mean spatial extension. 
All quantities are in expressed reduced units. 



At later times the dynamics will not further affect energy 
and fragment partitions, meaning that the comparison we 
show will be pertinent also for later times. 



Energy sharing and effective temperature 

The time dependent sharing of the internal energy 
Eint = Kint + V between the kinetic and interaction 
component is shown in Figure 2, which displays the av- 
erage disordered kinetic energy (Ki nt ) as a function of the 
internal energy for the different evolutions. The equilib- 
rium correlations for the initial and final state are also 
represented. 
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FIG. 2: (Color online) Average kinetic energy per particle as 
a function of the internal energy per particle of a confined 
LJ system in a harmonic trap with uj — O.Olto" 1 (full line) 
and uj — O.Uq 1 (dashed line) compared to LJ systems in free 
expansion at different total energies at times between 0to and 
tf = 250to (symbols). Emphasized symbols correspond to 
the time at which the confined and the free expanding system 
have the same mean spatial extension. 



We can see that the equilibrium sharing is well verified 
at the freeze-out time, except a slight underestimation 
at the lower energies. This result implies that the ki- 
netic energy at freeze-out, once the collective component 
is subtracted, can be used as a thermometer measuring 
a physically well defined temperature for the expanding 
system. 

The deviation at low energy is also interesting. At 
E/N — — le the system evolution is essentially a sim- 
ple evaporation of monomers and light fragments at the 
surface of a single excited condensed drop. The poor 
reproduction by the statistical ansatz (Q} of the drop ki- 
netic energy implies that the evaporation rate can only 
approximately be described in terms of an effective pres- 
sure , due to the irreducible time dependence of the 
evaporation process. 
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Statistical treatment of the expansion dynamics 

At first sight it may be surprising that the best ade- 
quacy to the equilibrium picture is obtained at the two 
highest energies, where the flow contribution is the most 
important, and the dynamics is the fastest (see Figure 
1). This counter- intuitive result can be understood if we 
consider that the structure (Q} for the microstate distri- 
bution is exactly preserved by the time evolution in the 
case of non-interacting particles or local interactions, as 
we now show 
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|23(. For a Hamiltonian system, the 
time dependence of any mean observation (A) is given 
by dt(A) = —({H,A}). If the system is non-interacting 
(H = K) this immediately gives 



(7) 



It is easy to show that these same relations hold for a 
non-Hamiltonian dynamics in the presence of a Boltz- 
mann collision integral accounting for local two body in- 
teractions. Let us consider an initial condition given by 
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FIG. 3: (Color online) Time dependence of the average 
square radius, radial momentum, and square momentum (all 
in reduced units) for the freely expanding system at four 
different energies. Symbols: numerical simulations, where the 
internal square momentum of the largest cluster recognized 
at asymptotic times is subtracted. Lines: free gas evolution 
eqs. (|15ll6H7p corresponding to the statistical ansatz eq. (|14p 
with an initial flow. 



the equilibrium distribution eq.(|3]) (with vlj = 0), im- 
posing an average kinetic energy and mean square radius 
through the Lagrange multipliers (3 and A = f3mu) 2 /2 
constraining respectively the observables (p 2 ) and (r 2 ). 
At any successive time the exact evolution of these ob- 
servables will be given by eqs.JT]). The minimum biased 
distribution fulfilling these time dependent constraints 
and additionally conserving the total energy, is given at 
any time t by 



P (n) (t) 



61E-H™ 
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■ exp 
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(8) 



(9) 



Using the associated dynamical equation (Liouvillc equa- 
tion for the ideal gas or Boltzmann equation for the col- 
lision problem), it is possible to show [2l[ that indeed 
eq. (JSj) is the exact evolution of eq. (jlj under the action of 
the ideal gas Hamiltonian. 

Eq. (J8J) can be interpreted as a radially expanding ideal 
gas in local equilibrium with a time dependent tempera- 
ture and pressure. Indeed the distribution can be written 
as 



P [n) (t) = 



exp 



S [E-H 



L.T 



W x (E,t) 



(10) 



N 



-£(*)£ 



2m 



where h = l/t represents a Hubblian factor. 

In the absence of an attractive inter-particle interac- 
tion the whole momentum distribution participates to 
the flow dynamics. As noticed above, in the LJ case a 
part of the initial kinetic energy is converted into internal 
energy of the clusters. Because of the spherical symme- 
try of the problem (%2i=iPi) = ^(Pr)> the two contribu- 
tions are decoupled N(p 2 )/2m = (K int ) + (Efi ow ) (see 
eas. (|5|6]> ). and we can expect eqs.((7|) to be still approxi- 
mately fulfilled, provided the internal contribution to the 
kinetic energy is subtracted. The time evolution of the 
average square radius, collective momentum, and square 
momentum are represented for different total energies in 
Figure 3. At the beginning of the evolution the strong 
interactions acting in these dense supercritical systems 
modify the dynamics with respect to the ideal gas or di- 
luted Boltzmann ansatz. However, starting from time 
t w 20i we can see that the trend predicted by eqs.([7]) 
is well verified. At t ~ 20io the system is still dense and 
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homogeneous in first approximation [l5| (see fig. 1) . We 
can then make the assumption that at t « 20to = t\ 
the initial condition (Q} still describes the observed dis- 
tribution with the extra constraint of a collective flow 
(p ■ r) (h) £ 0, 



cqs. (|15ll6H7[) shows that, even in the diabatic limit, as 
far as average global observables are concerned, the ex- 
pansion can be represented at any time as a statistical 
equilibrium in the local rest frame eq.©. 
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exp 
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i=l 



In the ideal or Boltzmann gas limit, the successive evo- 
lution of the distribution (fTTj) is given again by eq.©, 
with a modification of the time dependent constraints to 
account for the initial flow. For the purpose of getting 
analytical results it is easier to consider an initial condi- 
tion in the canonical ensemble: 



Pi±(ti) = y i -exp[-EA(p1 n) ) 2 -A 1 (rf ) ) 2 + 

(12) 
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Then the time dependent partition sum can be factorized 
Zs c ~ = z 1 ^ _ with 



1 



d 3 1 



d 3 P exp(-/3i (t)^--Ai (t)r- 



2m 



v\ {t)p- r 

2y / 7T77l 



(13) 



h 3 



(20 



1A1 



v\ 2 m 



-3/2 



The time dependent Lagrange parameters are given by 
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and the time dependent equations of state give for any 
Af = t — 1\ > the predicted evolution of the average 
observables: 
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d log zs <■ - 



6mAi 



r- = -- 



0/Mt) 

»Ai(t) 



2/3iAi — m^i 2 

6Ai 



(15) 



3/3i - GvtAt 



At 2 



2Pi\i — mv\ 2 



-{16) 



9 log 2^ - 3mui-6AiAt 
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The good adequacy displayed in Figure 3 between 
the numerical evolution and the ideal gas dynamics 



Deviations from equilibrium 

The good reproduction of time dependent average ob- 
servables by the equilibrium ansatz does not mean that 
the configurations explored by the expanding system co- 
incide with the equilibrium configurations. To look for 
deviations from equilibrium, we first show in Fig. 4 the be- 
havior of potential energy fluctuations. In the low energy 
liquid regime, as well as close to the liquid-gas transition, 
the freeze-out fluctuations are very close to the prediction 
of the equilibrium model. For the intermediate energies 
(E/N = Oe, E/N = le) however, these fluctuations are 
underestimated by the equilibrium calculation. 

This underestimation can be partly due to the energy 
conservation constraint. Indeed in the expansion of an 
isolated system the total energy is a constant of motion. 
This means that the potential energy can be converted 
both into internal kinetic energy, and into collective flow, 
inducing event-by-event fluctuations in both components. 
This effect is not considered in eq.([T]) where the flow is 
absent, and the conservation law therefore applies on the 
total internal energy. In other words, the equilibrium 
ansatz ^ misses flow energy fluctuations. 

To quantify this effect and explore if it can explain the 
fluctuation underestimation, we have plotted in Figure 
4 also the prediction of eq. fTTj) . Three models are then 
compared in this figure: the free expansion at time tf, 
and the equilibrium model in the local rest frame includ- 
ing (eq.(QT])) or not (eq.([T])) flow fluctuations. The three 
considered models are tuned to have in average the same 
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FIG. 4: (Color online) Potential energy per particle fluctua- 
tions for the free expanding, constrained without flow (eq.|T])), 
and constrained with flow (eq. (|lip ') systems. The abscissa 
gives the total energy per particle of the free expanding sys- 
tem in all cases. 
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internal energy (Ei nt ) and spatial extension (R 2 ). This 
internal energy is fixed in each event (n), K^ n > + V^ n > — 
Ei n t — constant, in the case of eq.([T]); while it can fluc- 
tuate in the free expansion and also for eq. flT]) , where 
KM + KM = 41 + E$ ow = E tot = constant. We 
can see from Figure 4 that indeed flow fluctuations can 
enhance the potential energy variance, but they cannot 
explain the deviation as they do not exhibit the correct 
energy dependence. Indeed at the highest energy eq. (|ll|) 
strongly overshoots the expanding simulation. 

To understand the origin of this deviation, we plot in 
Figure 5 the matter density and radial velocity profiles 
associated to the four different energies and the three 
models shown in Figure 4. The statistical simulations al- 
ways present a well pronounced density peak close to the 
center of mass of the system, corresponding to a single 
drop of decreasing size. The rest of the matter concen- 
trates on a second density peak around r rs 15ct, where 
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FIG. 5: (Color online) Upper part: average number of parti- 
cles located in a spherical shell of radius r as a function of r, 
calculated from the center of mass of the system. Lower part: 
average velocity in the radial direction in the same spherical 
shells as the upper part. The 4 energies ( — 1.0,0.0, 1.0 and 
2.0e) and the three models (free expansion, equilibrium in 
the local rest frame with and without flow fluctuations) are 
displayed. 
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FIG. 6: (Color online) Distributions of potential energy (left 
panel) and of the size of the biggest MST fragment (right 
panel) for the highest energy considered (2.0e) and the three 
models (free expansion, equilibrium in the local rest frame 
with and without flow fluctuations). 



most of the clusters are located. 

The spatial distribution of matter for the freely ex- 
panding system is close to this picture at the lowest en- 
ergy under analysis. However the presence of a heavy 
drop close to the center of mass modifies the shape of the 
radial flow, which is very far away from the self-similar as- 
sumption of eq. (TTTj) . The Hubblian flow in eq. (fTTj) comes 
from the ideal-gas equations of motion eq.([7]). We can see 
from Figure 5 that the inter-particle LJ interaction can- 
not be neglected in the dense part of the system, which 
leads to an important deviation from self-similarity We 
expect that a better description of the expansion dynam- 
ics would be obtained considering for the flow dynamics 
in ea. fTTj) clusters instead of particles degrees of freedom, 
as in the standard Fisher picture of condensation 24| . 

At higher energies, the distribution of the free expand- 
ing system is completely different from the equilibrium 
picture. The central drop progressively disappears lead- 
ing to a smoother matter distribution along the radial 
direction. From previous studies [l6| . we know that the 
resulting large bump corresponds to a high multiplicity of 
approximately equal sized clusters. The presence of clus- 
terized matter leads to a deviation from the Hubblian 
flow even at the highest energy. 

More information on the specific configurations ac- 
cessed by the flow dynamics can be obtained from Figure 
6. This figure displays the distribution of the potential 
energy and of the size of the l arge st fragment recognized 
through the MST algorithm [16(. It confirms that the 
freely expanding system explores configurations which 
are very different from the equilibrium case. As already 
observed in Figure 2, the global energy sharing between 
potential and kinetic energy of the free expanding sys- 
tem is consistent with equilibrium, and consequently the 
potential energy distributions are close. However the as- 
sociate partitions are very different and tend to be more 
fragmented in the free expanding shown by the 

fact that the largest cluster size distribution is broader 
and peaked at a lower value. At the energy shown in the 
figure, the system is close to the liquid-gas phase transi- 
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tion. If flow fluctuations are allowed, as in eq.(fTTj). some 
gas partitions can then be explored, leading to the bi- 
modal distributions observed in Figure 6. The fact that 
a statistical treatment of flow can lead to bimodality was 
recently observed in ref . [23| . The presence of bimodal dis- 
tributions explain the severe overestimation of the free 
expanding system fluctuations shown in Figure 4. In- 
deed, for these distributions to be accessed during the 
actual time evolution of the expansion a potential energy 
barrier has to be overcome, and this does not seem to 
occur easily for the free expanding system. As it was al- 
ready shown in ref. [25} , flow acts as a heat sink preventing 
the exploration of gas configurations. The largest cluster 
distribution shows that most partitions of the freely ex- 
panding system present an intermediate degree of frag- 
mentation between liquid and gas. Such partitions are 
metastable or unstable at equilibrium, but can be ac- 
cessed in the free expansion due to the short time scale 
of the dynamics. 

CONCLUSIONS 

To conclude, in this paper we have compared the dia- 
batic expansion dynamics of a Lennard- Jones system, ini- 
tially confined in a harmonic oscillator and subsequently 
expanding freely in the vacuum, with a statistical ansatz 
in the hypothesis of a purely Hubblian flow. This hy- 
pothesis is exact in the limiting case of a non-interacting 
system or a Boltzmann dynamics. For our strongly inter- 
acting system, the presence of finite range two body in- 
teractions is known (2l| to modify the Hubblian approxi- 
mation introducing non-self-similar flow components. In 
a future work, it will be very interesting to explore the 
adequacy of a more sophisticated statistical ansatz in- 
cluding non-self-similar flows, to reproduce the dynamics 
of the expansion. In the present paper, we have lim- 
ited ourselves to the self-similar approximation which re- 
veals to be accurate enough to reasonably describe the 
mean value of global one-body observables at all times. 
As soon as more sophisticated observables are exam- 
ined, discrepancies arise. In the self-similar statistical 
ansatz, flow does not modify the partitions in configura- 
tion space, but it acts as a heat bath allowing important 
energy fluctuations and the exploration of the unbound 
gas phase. In the diabatic dynamics, such configurations 
are never reached and metastable highly clusterized par- 
titions dominate. It is interesting to observe that quali- 
tatively similar behaviors have recently been observed in 
an analysis of nuclear multifragmentation data [26j], by 
means of a detailed comparison between the fragmenta- 
tion of central and peripheral collisions. 
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